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Abstract 

Background: Meat quality depends on skeletal muscle structure and metabolic properties. While most studies carried on 
pigs focus on the Longissimus muscle (LM) for fresh meat consumption, Semimembranosus (SM) is also of interest because of 
its importance for cooked ham production. Even if both muscles are classified as glycolytic muscles, they exhibit dissimilar 
myofiber composition and metabolic characteristics. The comparison of LM and SM transcriptome profiles undertaken in 
this study may thus clarify the biological events underlying their phenotypic differences which might influence several meat 
quality traits. 

Methodology/Principal Findings: Muscu\ar transcriptome analyses were performed using a custom pig muscle microarray: 
the 1 5 K Genmascqchip. A total of 3823 genes were differentially expressed between the two muscles (Benjamini-Hochberg 
adjusted P value <0.05), out of which 1690 and 2133 were overrepresented in LM and SM respectively. The microarray data 
were validated using the expression level of seven differentially expressed genes quantified by real-time RT-PCR. A set of 
1047 differentially expressed genes with a muscle fold change ratio above 1.5 was used for functional characterization. 
Functional annotation emphasized five main clusters associated to transcriptome muscle differences. These five clusters 
were related to energy metabolism, cell cycle, gene expression, anatomical structure development and signal transduction/ 
immune response. 

Conclusions/Significance:lV\\s study revealed strong transcriptome differences between LM and SM. These results suggest 
that skeletal muscle discrepancies might arise essentially from different post-natal myogenic activities. 
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introduction 

Pork is one of the most widely eaten meats in the world. 
Breeding programs aiming at improving pig production efficiency 
through increased growth rate and lean meat content and 
decreased fatness have also affected some meat quality traits 
playing an important role in consumer acceptance of pork like 
water holding capacity, color, intramuscular fat (IMF) content and 
tenderness [1]. Meat quality is a complex trait which depends on 
the interactive effects of pig genotype, environmental conditions, 
pre-slaughter handling and slaughtering procedure [2]. The 
skeletal muscle structure and metabolic characteristics which 
determine cellular and molecular events occurring during muscle 
to meat transformation are of the utmost importance for meat 
quality determination. Skeletal muscle is a heterogeneous tissue 
composed of myofibers, adipose, connective, vascular and nervous 
tissues. Myofibers differ by their molecular, structural, contractile 



and metabolic properties according to which they are classified as 
slow-twitch oxidative (type I), fast-twitch oxido-glycolytic (type IIA) 
and fast-twitch glycolytic (type IIB). Red or white muscles are also 
determined according to their fiber type composition. Red muscles 
are composed of high percentage of slow-twitch oxidative fibers 
whereas white muscles contain a major proportion of fast-twitch 
glycolytic fibers [3]. Longissimus and Semimembranosus - two white 
skeletal muscles - are consumed in different forms: fresh for LM 
(loin) or after processing for SM (ham). Both muscles are classified 
as glycolytic even if shght differences have been described in their 
myofiber composition (higher proportion of type Ila myofiber and 
lower proportion of Type lib myofiber in SM) and metabolic 
properties (higher oxidative capacity in SM) [4—7] . Transcriptome 
analysis might be useful to identify transcriptional signatures 
associated with meat quality traits which could thus be selected as 
biomarkers in selection programs [8-12]. However, pig transcrip- 
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tome studies are mainly focused on LM even if gene expression 
variability between muscles could affect muscle development, meat 
quality and hence the choice of meat processing [13]. 

The aim of this study was to better characterize LM and SM 
gene expression profiles in order to investigate the biological 
events underlying their distinct metabolic and contractile proper- 
ties. 

Results 

Comparison of Gene Expression Profiles between 
Longissimus and Semimembranosus Muscles 

Gene expression microarray analysis was conducted on 180 
muscle samples (90 LM, 90 SM). Comparison of LM and SM 
muscle transcriptome was achieved using the "GenmascqChip", a 
15 k pig skeletal muscle microarray. Raw data sets were clu-ckcd 
for quality criteria. The 10753 remaining probes were considered 
significandy expressed in both muscles. We observed a strong 
muscle effect on gene expression with 5582 (37%) of probes being 
differentially expressed between LM and SM (adjusted P value S 
0.05). As shown in Figure 1, fold change (FC) ratios varied from 
1.1 to 15 and were for the most part quite low with median values 
<1.5 in the two muscles. These 5582 difiTerentially expressed 
probes corresponded to 3823 annotated genes, with 1690 and 
2133 genes overrepresented in LM (Table SI) and SM (Table S2), 
respectively. A set of 2402 difierentiaUy expressed probes (1603 
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Figure 1. Gene expression ratio between muscles. Muscle fold 
change ratio is expressed as the expression ratio of Longissimus (LIVl) to 
Semimembranosus (SM) samples when genes are highly expressed in 
Longissimus and as the expression ratio of SM to LM samples when 
genes are highly expressed in Semimembranosus. 
doi:1 0.1 371/journal.pone.0096491 .gOOl 



annotated genes) with a muscle FC ratio above 1.5 was considered 
as biologically relevant and selected for functional analysis. 

The ten most differentially expressed and informative genes [i.e. 
with at least one associated gene ontology' (GO) biological process 
(BP) term] are shown in Table 1 for LM (5.3<FC< 15.3) and in 
Table 2 for SM (3.4<FCS8.1). Among the ten genes strongly 
expressed in the LM, three are involved in gene expression: 
poly(rC) binding protein 2 (PCBP2), microphthalmia-associated 
transcription factor [MITF) and zinc finger and BTB domain- 
containing protein 16 [ZBTB16). Three other genes are involved in 
metabolism: ATP synthase, H+ transporting, mitochondrial Fo 
complex, subunit E {ATPSTj, ADAM metallopeptidase with 
thrombospondin type 1 motif, 8 {ADAMTS8) and kinase D- 
interacting substrate, 220 kDa {nDINS220). ADAMTS8 and 
ZBTB16 are also related to negative regulation of cell prolifera- 
tion. Two genes are involved in muscle development: interferon- 
related developmental regulator 1 (IFRDl) and ryanodine receptor 
1 [RYRl) which is also involved in muscle contraction and calcium 
ion transport. Last, one gene is related to cell-cell signaling: nudix 
(nucleoside diphosphate linked moiety X)-type motif 3 (NUDT3) 
and one gene in DNA replication and DNA repair: REVS-Uke, 
polymerase (DNA directed), zeta, catalytic subunit {REV3L). For 
the SM, four genes are involved in gene c'xpression: spleen focus 
forming virus proviral integration oncogene (SPIT), nuclear 
receptor subfamily 2, group C, member 2 {NR2C2), histone cluster 
1, H2ab (HIST1H2AB) and tenascin C (rjVC) which is also related 
to positive regulation of cell proliferation. Two genes are involved 
in muscle contraction: protein phosphatase 1, regulatory subunit 
12B {PPPlR12Bj and myosin, heavy chain 1 1 {MTHll). Last, one 
gene is involved in wounding and inflammation (acid phosphatase 
5, tartrate resistant, ACP3), one gene in water transport (aquaporin 
4, AQP4), one gene in cellular component movement (kinesin 
family member C2, KIFC2) and one gene in cell adhesion (secreted 
phosphoprotein 1, SPPl). 

Quantitative RT-PCR Validation of Microarray Analysis 

Seven target genes, including four genes overrepresented in LM 
{ADAMTS8, ALDOA, CPTIB and RYRI) and 3 genes overrepre- 
sented in SM [CEBPA, DGAT2 and TGFBI) were analyzed by real 
time RT-qPCR. These genes were selected to represent the 
variation of FC ratio observed across the set of 1603 differentially 
expressed genes with a muscle FC ratio above 1.5. As shown in 
Figure 2, the comparison of FC ratios between microarray and 
RT-qPCR technologies provided similar FC direction between 
these two methods. However, FC values were much less consistent 
between methodologies for gene overexpressed in LM than in SM 
samples. 

Functional Analysis 

To identify the biological events to which differentially 
expressed genes product contributes, we used GO BP annotation. 
A set of 2402 differentially expressed probes with a muscle FC 
ratio above 1.5 was selected for functional analysis. They 
corresponded to 1603 human orthologous genes. Among them, 
1047 were associated with at least one GO BP term and were 
clustered according to their semantic similarities using these terms 
(Figure 3). Cluster compositions are shown in Table S3. 

Five clusters related to energy metabolism (cluster 1 including 
142 genes), cell cycle (cluster 2, 175 genes), gene expression (cluster 
3, 127 genes), anatomical structure development (clusters 4, 480 
genes) and cell communication/immune response (cluster 5, 123 
genes) were identified. For each cluster, some relevant GO BP 
terms and pathways (KEGG and WikiPathways) are presented in 
Table 3 and Table 4. Full details of enriched biological processes 
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Table 1. Genes overexpressed in Longissimus muscle (n = 10). 



Symbol^ 



Cluster^ 



FC^ 



P-value 



Associated GO BP terms 



ADAMTS8 



RYR1 



REV3L 



rPRDI 



PCBP2 



NUDT3 5 
ATP5I 1 



KIDINS220 



MITF 



ZBTB16 



13.6 



8.3 



6.4 



6.3 
5.8 



5,3 



5,3 



<1E" 



<1E" 



<1E" 



<1E" 



<1E" 



<1E" 
<1E" 



<1E" 



<1E" 



<1E" 



GO:0008285"-Negative regulation of cell proliferation 
GO:0006508~Proteolysis 
GO:000681 6— Calcium ion transport 
GO:0048741 —Skeletal muscle fiber development 
GO:0006936~Mu5de contraction 
GO:0006261~DNA-dependent DNA replication 
GO:0006281~DNA repair 
GO:0007518~Myoblast cell fate determination 
GO:0042692~Mu5cle cell differentiation 
GO:0007527— Adult somatic muscle development 
GO:0008380~RNA splicing 
GO:001 0467— Gene expression 
GO:0016071~mRNA metabolic process 
GO:0007267~Cell-cell signaling 
GO:0022904— Respiratory electron transport chain 
GO:0006200~ATP catabolic process 

GO:0042776— Mitochondrial ATP synthesis coupled proton transport 

GO:0000186~ Activation of MAPKK activity 

GO:004801 1 —Nerve growth factor receptor signaling pathway 

GO:0007275— Multicellular organismal development 

GO:0045893~Positive regulation of transcription, DNA-dependent 

GO:0006915~Apoptotic process 

GO:0008285— Negative regulation of cell proliferation 

GO:0045893~Positive regulation of transcription, DNA-dependent 

GO:0045892— Negative regulation of transcription, DNA-dependent 



Only genes with at least one associated GO BP term are presented in the table. 

^Differentially expressed genes were clustered using GO BP terms semantic similarity between genes as distance, to group functionally similar genes together. 
^Fold Change is expressed as the expression ratio of Longissimus to Semimembranosus samples. 
''Benjamini and Hochberg adjusted P value. 

^Unique identifier and gene ontology term in the GO database (http://www.geneontology.org/). 
doi:l 0.1 371 /journal.pone.0096491 .tOOl 



and pathways, enrichment score, adjusted P-value and number of 
gene present in each cluster are reported in the Table S4 and 
Table S5. Cluster 1 comprised 98 genes highly expressed in SM 
and 44 in LM. Significantly enriched GO BP terms (P-value < 
7.4E~'", enrichment score (ES): 1.4 to 15.2) and pathways (P- 
value <1.8E— 02, ES: 3.2 to 38) were mainly related to energy 
metabolism. Cluster 1 genes were assigned to several enriched 
biochemical pathways including "Electron Transport Chain", 
"Oxidative phosphorylation", "Glycolysis and Gluconeogenesis", 
"Fatty Acid Beta Oxidation" and "Citrate cycle (TCA cycle)". 
The GO BP term "generation of precursor metabolites and 
energy" with the highest P-value, was associated with 39 genes 
encoding five mitochondrial electron transfer chain complex 
subunits that were mainly expressed in SM. Succinate dehydro- 
genase complex, subunit A, flavoprotein (Fp) [SDHA) and genes of 
long chain fatty acid metabolism (acyl-CoA dehydrogenase, very 
long chain, ACADVL) were also more expressed in SM whereas 
solute carrier family 25, member 27 {SLC25A27 also known as 
uncoupling proteins 4 UCP4) and solute carrier family 25, member 
14 [SLC25A14 also known as uncoupling proteins UCPSj, 
phosphorylase kinase, alphal and beta {PHKAl, PHKB) and 
phosphoenolpyruvate carboxykinase 1 [PCKT) were overexpressed 
in LM. Cluster 2 included 73 and 102 highly expressed genes in 



SM and LM, respectively. Enriched GO BP terms (P-value < 
e.eE"""^, ES: 1.1 to 7.5) and pathways (P-value <1.6E-02, ES: 
2.7 to 7.6) were related to cell cycle process. "Cell cycle" and 
"Ubiquitin mediated proteolysis" were the most important 
enriched pathways associated with cluster 2. Genes overrepre- 
sented in SM were mainly linked to Gl phase: cyclin D2 and D3 
{CCND2, CCND3). Genes overrepresented in LM were related to 
the control of cell cycle checkpoint, GO/Gl, Gl/S, S/G2 and 
G2/M transition and M phase: anaphase promoting complex 
subunit 1 and 4 {ANAPCl, AMAPC4), cell division cycle 26 and 27 
(CDC26, CDC27) as well as DNA rephcation and DNA repair 
process. Cluster 3 contained 43 and 84 highly expressed genes in 
SM and LM, respectively. Cluster 3 enriched GO BP terms (P- 
value <1.6E ES: 1.5 to 11.7) were related to gene expression. 
Significantly enriched pathways were related to "mRNA process- 
ing" (P-value = 1.5E— 10, ES = 9.3) and "Spliceosome" (P-va- 
lue =1.7E- 12, ES=10.9). Li this cluster, two of the four 
myogenic regulatory factors, myogenic differentiation 1 [MYODT) 
and myogenic factor 6 (MTF6 also known as MRF4) were 
overrepresented in SM. Cluster 4 was the biggest one with 288 
and 192 genes overexpressed in SM and LM, respectively. 
Enriched GO BP terms were mainly related to anatomical 
structure development (P-value <1.7E ES: 1.5 to 4.4). Among 
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Table 2. Genes overexpressed in Semimembranosus nnuscle (n = 10). 



Symbol^ 



Cluster^ FC' 



P-value 



Associated GO BP terms 



SPI1 



NR2C2 



PPP1R12B 



TNC 



SPP1 



HIST1H2AB 
AQP4 



ACP5 



KIFC2 
MYHll 



5.6 



53 



4.3 

3.9 



3.7 



35 
3,4 



<1E" 



<1E" 



<1E" 



<1E" 



<1E" 



<1E" 
<1E" 



<1E" 



<1E" 
<1E" 



GO:0045893~Po5itive regulation of transcription, DNA-dependent 

GO:0045892— Negative regulation of transcription, DNA-dependent 

GO:0045814~Negative regulation of gene expression, epigenetic 

GO:0006355~Regulation of transcription. DNA-dependent 

GO:0030154~Cell differentiation 

GO:0010467"'Gene expression 

GO:0006937'-Regulation of muscle contraction 

GO:0007165~Signal transduction 

GO:0008284~ Positive regulation of cell proliferation 

GO:0007528"'Neuromuscular junction development 

GO:0010628'-Positive regulation of gene expression 

GO:0007155~Cell adhesion 

GOiOOOl 649~Osteoblast differentiation 

GO:0001 503~Ossification 

GO:0006334~Nucleosome assembly 

GO:000681 0~Transport 

GO:0006833~Water transport 

GO:0050891 —Multicellular organismal water homeostasis 
GO:0060349"-Bone morphogenesis 

GO:0050728— Negative regulation of inflammatory response 

GO:000701 8~Microtubule-based movement 

GO:0030241 —Skeletal muscle myosin thick filament assembly 

GO:0048251 -Elastic fiber assembly 

GO:0006936-Muscle contraction 



^Only genes with at least one associated GO BP term are presented in the table. 

^Differentially expressed genes were clustered using GO BP terms semantic similarity between genes as distance, to group functionally similar genes together. 
^Fold Change is expressed as the expression ratio of Semimembranosus to Longissimus samples. 
''Benjamini and Hochberg adjusted P value. 

^Unique identifier and gene ontology term in the GO database (http://www.geneontology.org/}. 
doi:l 0.1 371 /journal.pone.0096491 .t002 



enriched pathways found associated with cluster 4, "ECM- 
receptor interaction" (P-value = 7. IE — 07, ES = 5.9) and "Focal 
adhesion" (P-value = LIE — 6, ES = 3.4) had the highest P-value. 
Genes overrepresented in LM were implicated in cell division, 
chromosomal organization, structural maintenance of sarcomere 
and sarcoplasmic protein: nebulin (MEB), titin (TZN), RYRl and 
triadin [TRDM). Genes overrepresented in SM were involved in 
cell migration, cell surface and extracellular matrix (ECM): 
cadherin 2, type 1, N-cadherin {CDH2), CD44 molecule {CD44), 
caveolin 3 {CAV3) and code for major constituent of the contractile 
apparatus: myosin, heavy chain 3, 8, 9 and 11 {MYH3, MYH8, 
Mm9 and MYHll), troponin I type 3 {TNNB) and troponin T 
type 2 {TMNT2). Differentially expressed genes that might control 
muscle size were either more expressed in LM, foUistatin {FST) and 
myostatin (MSTN) or in SM, insulin-like growth factor 1 
(somatomedin C) {IGFl), transforming growth factor, beta 1 and 
3 (TGFBl and TGFB3). Last, cluster 5 contained 90 overexpressed 
genes in SM and 33 overexpressed genes in LM. Significantly 
enriched GO BP terms (P-value <2E-05, ES: 1.7 to 16.3) and 
biological pathways (P-value <1.E — 12, enrichment score (ES): 2.6 
to 39) were related to cell communication and inflammatory 
immune response. Interestingly, this cluster highlighted SM 
overexpressed genes also involved in the regulation of develop- 
mental and myogenesis signaling pathways or in the muscle 



regeneration process: fibroblast growth factor 18 [FGFIS), 
members of the Notch signaling pathway [N0TCH3 and delta- 
like 4, DLL4), chemokine (C-C motif) hgand 21 and 23 {CCL21 
and CCL23) and complement component and factor (complement 
component 1, q subcomponent, C chain, CIQC; complement 
component 1 , r subcomponent, CIR; complement component 1 , s 
subcomponent, CIS; complement component 3, C3; complement 
component 4A, C4A and complement factor B, CFB). 

Discussion 

Our objective was to clarify the biological events which could 
explain the muscle phenotypic differences reported in the 
literature between the LM and SM [3,4,14]. Since skeletal muscle 
is a heterogeneous tissue, transcriptome analysis of skeletal muscle 
may reflect mRNA composition of various cell types existing in 
this tissue. However, we assumed that myofibers are the main 
skeletal muscle component and that comparison between muscles 
is informative. The custom GenmascqChip [15] used in this study 
allows the analysis of 10753 probes and the identification of 5582 
diflFerentially expressed probes between LM and SM demonstrat- 
ing that the GenmascqChip is a powerful tool to study pig muscle 
gene expression in order to gain a better understanding of muscle 
physiology. Furthermore, directions of differential gene expression 
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Figure 2. Validation of seven microarray differentially expressed genes between Longissimus (LM) and Semimembranosus (SM) 
muscles by quantitative RT-PCR. mRNA level is expressed using arbitrary units. Quantitative RT-PCR expression levels (LM = 8, SIVI = 8) were 
normalized to the expression of beta 2 microglobulin {B2M), TATAbox binding protein (7BP) and 18S using geNorm algorithm. Microarray adjusted 
means for LM and SM (LM = 90; SM = 90) were calculated using least square means for the muscle effect. Data are expressed as means±s.d. Statistical 
significances are reported below the plot as Benjamini and Hochberg adjusted P-value for microarray data and as Student t-test P value for q RT-PCR. 
Fold change ratio is expressed as the expression ratio of LM to SM when genes are overrepresented in LM and as the expression ratio of SM to LM 
when genes are overrepresented in SM. 
doi:10.1371/journal.pone.0096491.g002 



FC observed in microarray analysis have been validated by 
quantitative PGR analyses using seven differentially expressed 
genes. We have considered that FC direction agreement between 
the two methods validates the differentially expressed genes set. 
Besides, there are few studies comparing gene expression of 
contrasted skeletal muscles in pigs, most of them focus on LM and 
none of them compared LM and SM [16,17]. The number of 
genes found differentially expressed between these two muscles is 
surprisingly high. In fact, Hornshoj et al. [17] described a similar 
expression pattern between these two muscles, while the compar- 
ison of red and white skeletal muscles that are much more 
contrasted than LM and SM led to far fewer differentially 
expressed genes in pigs [16,18] or mice [19]. However, different 
experimental conditions such as microarray platform technology, 
whole genome vs. focused microarray, sample size, samples 
pooling, FC threshold might account for this discrepancy. Several 
studies have compared gene expression level from different 
microarray technologies and relate divergence across the data 
generated [20-23]. Stretch et al. [24] have studied the effect of the 
sample size on differentially expressed gene discovery. They 
studied muscle gene expression on 134 samples (69 males, 65 
females) and found that using sample of n = 1 0 (5 males, 5 females) 
results in no significant genes at P-values <0.0001, whereas larger 
sample size n= 120 (60 males, 60 females) identifies 472 
differentially expressed genes at the same P-value cutoff. Anyway, 
this new finding reinforces the importance of gene expression 
variability between muscles which could affect muscle develop- 
ment and hence meat quality [13]. 

Microarray experiments result in list of hundred to thousand 
differentially expressed genes and the main objective of functional 



data analysis is to determine relevant biological interpretations. In 
this context, hierarchical clustering is often performed using gene 
expression correlation coefficient matrix as distance considering 
that co-expressed genes share the same biological processes. 
Biological knowledge is then used to identify enriched biological 
processes in each gene cluster [25]. Using this approach, we 
obtained two large clusters corresponding to over- and underrep- 
resented genes which led to dozens of dissimilar enriched terms 
(data not shown). Furthermore, biological pathways are mosdy 
controlled by the balance between up and downregulations. 
Performing functional analysis separately for up and downregu- 
lated genes list might result in loss of biological information since 
genes involved in the same pathway could have been assigned in 
different set. This partial information may leads to misinterpre- 
tation of differentially expressed gene list. Some pathways might 
then have been discarded because their enrichment value was 
deemed insufficient, whereas gene involved in the regulation of 
this pathways were present in up and downregulated genes list. To 
avoid this and create meaningful clusters, we used semantic 
similarity of GO BP terms to group functionally similar genes 
together. Wang's metrics [26] was chosen over the information 
content-based semantic similarity measures because the latter 
require a rehable corpus in order to compute GO terms 
frequencies, and such a corpus does not exist for moderately 
studied species such as Sus scrofa. We successfully identified five 
functional clusters including both over and underrepresented 
genes. This approach was well suited to the size of our data set 
(around 1600 genes). However, the hierarchical clustering 
algorithm led to exclusive classification and we assume that 
clusters cannot overlap whereas genes may be involved in several 
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Figure 3. Hierarchical clustering of differentially expressed genes according to their GO BP terms semantic similarity. Annotated 
differentially expressed genes with a muscle fold change above 1.5 were clustered based on their functional annotation (GO BP) semantic similarity. 
Hierarchical clustering was performed using "1 -semantic similarity" as distance between two genes (similar genes have a distance close to zero) to 
identify clusters of genes sharing BP terms. Five clusters were identified. Cluster 1 comprised 98 genes highly expressed in SM and 44 in LM. Cluster 2 
included 73 highly expressed genes in SM and 102 in LM. Cluster 3 contained 43 highly expressed genes in SM and 84 in LM. Cluster 4 comprised 288 
genes overexpressed in SM and 192 in LM. Cluster 5 involved 90 overexpressed genes in SM and 33 overexpressed genes in LM. 
doi:1 0.1 371/journal.pone.0096491 .g003 



biological processes. Thus most but not all relevant GO BP terms 
and pathways were highlighted with this procedure. The five main 
relevant biological networks associated to skeletal muscle differ- 
ences were "energ\' metabolism", "cell cycle", "gene expression", 
"anatomical structure development" and "cell communication/ 
immune response". Some examples of differentially expressed 
genes wiU be discussed in relation to energy metabolism and 
myogenic progenitor cells recruitment and sarcomerogenesis 
which composed steps of myogenesis process leading to the 
formation and growth of myofibers. The last part will discuss 
contrasted results in relation to muscle regeneration process. 

Energy Metabolism 

Our functional analysis identified "energy metabolism" as one 
of the most relevant biological pathway associated to LM and SM 
differentially expressed genes set. SM (n-erexpressed genes were 
related to mitochondrial fatty acid beta-oxidation pathway 
(ACSF3, ACADVL, ACADS and HADHA), citric acid cycle {AC02 
and SDHA) and the five mitochondrial respiratory chain complex: 
NADH dehydrogenase (ubiquinone) subunits {MT-MD3, MT-ND6, 
MDUFA3, JVDUFA8, MDUFA9, MDUFAIO, NDUFAU and 
NDUFVl), succinate dehydrogenase subunits (SDHA), ubiquinol- 
cytochrome c reductase complex subunits [CTC and UQCRCl), 



cytochrome c ()xidas(; subunits (00X412, C0X8A and MT-COl) 
and ATP synthase subunits {ATP5A1 and ATP5D). On the other 
hand, LM overexpressed genes were related to glycogenolysis 
regulation (PHKAl and PHKB), pyruvate metabolism pathways 
(PCK'l) and uncoupling protein (IJCP4 and UCP5). These results 
suggest on the one hand a higher mitochondrial oxidative activity 
in SM than in LM while on the other hand a limited usage of 
oxidative phosphorylation through uncoupling protein overex- 
pression and a predominant usage of the anaerobic glycolytic 
pathway in LM. These results are consistent with and refine 
previous knowledge on SM and LM metabolic characteristics. In 
fact, these two glycolytic muscles are predominantiy composed of 
fast-twitch t^pe II fibers and low level of slow-twitch type I fibers. 
However, SM is composed of highest percentage of intermediate 
fast-twitch type Ila myofibers and exhibited higher oxidative 
capacity than LM [3,4,6]. 

Myogenesis Process 

Although mature myofibers are postmitotic cells, functional 
enrichment analysis highlighted cell cycle, gene expression and 
muscle development as important features to characterize 
contrasted LM and SM expression profiles. SM overexpressed 
genes were related to satellite cells activation [IGFl, FGF18), cell 
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Table 3. Relevant biological processes significantly enriched in clustered differentially expressed genes. 





Cluster^ 


ES^ 


Specific GO term^ 


nG" 


P-value' 


1 




Energy metabolism 


142 






7.7 


GO:0006091 —generation of precursor metabolites and energy 


39 


<1E-'^ 




4.7 


GO:0006629~lipid metabolic process 


45 






5.1 


GO:00551 14~oxidation reduction 


41 


<1E-" 




9.2 


GO:0015980~energy derivation by oxidation of organic compounds 


24 


<1E-'^ 




10.3 


GO:a045333~cellular respiration 


18 


<1E-'^ 


2 




Cell cycle 


175 






2.8 


GO:0044267~cellular protein metabolic process 


110 


<1E-^^ 




2.5 


GO:0006464~ protein modification process 


54 


1.9E-™ 




4.4 


GO:0006281~DNA repair 


20 


3.7E-" 




6.8 


GO:0051439~regulation of ubiquitin-protein ligase activity during mitotic cell cycle 


9 


4.3E-°= 




2.3 


GO:0007049~cell cycle 


32 


47E-°= 


3 




Gene expression 


127 






5 


GO:0016070~RNA metabolic process 


113 


<1E-'^ 




3.4 


GO:0010467~gene expression 


119 


<1E-'^ 




7.7 


GO:0006396~RNA processing 


56 


<1E-' = 




10.2 


GO:0008380~RNA splicing 


42 


<1E-'^ 




3.1 


GO:0045449~regulation of transcription 


66 


<1E-'^ 




2.8 


GO:0010468~regulation of gene expression 


68 


<1E-'^ 


4 




Anatomical structure development 


480 






2.1 


GO:0009653~anatomical structure morphogenesis 


87 


3.4E-'"' 




2.5 


GO:0007155~cell adhesion 


60 


1.3E-'" 




1.8 


GO:0065008~regulation of biological quality 


97 


3.4E-'»' 




3.2 


GO:0003012~muscle system process 


27 


1.6E-" 




2.8 


GO:0007517~muscle organ development 


32 


1.9E-'* 


5 




Cell communication/immune response 


123 






3.1 


GO:0007154~cell communication 


96 


<1E-'^ 




3.2 


GO:0007165~signal transduction 


89 


<1E-'^ 




1.8 


GO:0050794~regulation of cellular process 


97 


<1E-'^ 




7.3 


GO:0006955~immune response 


39 


<1E-'^ 




6.2 


GO:0006954~ inflammatory response 


17 


1.9E-''= 



^Differentially expressed genes were clustered using GO BP terms semantic similarity between genes as distance, to group functionally similar genes together. 
^Cluster enrichment score (ES). 

^Unique identifier and gene ontology term in the GO database (http://www.geneontology.org/}. 
^Hq, number of genes in the category. 
^Benjamini and Hochberg adjusted P value. 
doi:l 0.1 371 /journal.pone.0096491 .t003 



cycle control at Gl phase [CCND2, CCND3 and cyclin-dependent 
kinase inhibitor, CDKNIB) and myoblast determination {MTODl, 
MRF4). On the other hand, LM overexpressed genes were 
involved in the negative regulation of satellite cells activation 
(MSTjV, FST), and ceU cycle progression through Gl/S, S/G2 or 
G2/M transition and in M phase. Hormonal control of satellite 
cells activation involved different growth factors including insulin- 
like growth factor I, which is a well-known hypertrophy factor 
acting on muscle mass and fibroblast growth factor [27,28]. In 
fact, insulin-Hke growth factor I induced myogenesis by activating 
satellite cells and promoting proliferation, differentiation and 
fusion with existing myoblast [27]. On the other hand, myostatin 
and its antagonist foUistatin, overexpressed in LM, are both 
involved in the main signaling pathway that negatively regulates 
satellite cells activation [29,30]. MSTN\s< expressed in satellite cells 
and act on cell cycle progression to maintain the Gl resting state 



(GO) and limit muscle growth by inhibiting satellite cells activation 
and proliferation [31]. FoUistatin antagonize myostatin inhibitory 
activity by direct protein interaction. Balance between foUistatin 
and myostatin limit the recruitment of sateUite ceUs [29]. 

Once activated, sateUite cells proliferate before undergoing 
myogenic differentiation. Cells proliferation relies on kinases or E3 
ubiquitin-protein Kgases that regulate activity or stability (ubiqui- 
tination and subsequent proteasomal degradation) of key ceU-cycle 
control proteins. Interestingly, we have identified "CeU cycle" and 
"Ubiquitin mediated proteolysis" among enriched pathways 
associated with cluster 2. Among E3 ubiquitin-protein ligases, 
Anaphase promoting complex/cyclosome (APC/C) genes 
lAMAPCl, ANAPC4, CDC26, CDC27) were overexpressed in LM. 
APC /C is a key regulator of the eukaryotic ceU cycle acting on 
GO/Gl transition, through S and G2 phase, and during mitoses to 
ensure proper and correct succession of ceU cycle key events [32- 
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Table 4. Relevant biological pathways significantly enriched in clustered differentially expressed genes 





Cluster^ 


ES^ 


Pathway name 


Pathways 


nG' 


P-value* 


1 




Energy metabolism 




142 






7.7 


Metabolic pathways 


KEGG 


84 


<1E-" 




15.7 


Electron Transport Chain 


Wiki pathways 


23 


<1E-" 




12.9 


Oxidative phosphorylation 


KEGG 


22 


<1E-" 




13.1 


Glycolysis and Gluconeogenesis 


Wiki pathways 


8 


1.7E-" 




7.5 


Fatty Acid Beta Oxidation 


Wikipathways 


8 


1.0E-°= 




13.2 


Citrate cycle (TCA cycle) 


KEGG 


5 


4.6E-''= 


2 




Cell cycle 




175 






6.6 


Cell cycle 


KEGG 


11 






5.3 


Ubiquitin mediated proteolysis 


KEGG 


11 


4.8E-°= 




6.4 


p53 signaling pathway 


KEGG 


6 


I.OE-"^ 




5.1 


RIbosome 


KEGG 


7 


1.1 E-o^ 




5.8 


Gap junction 


KEGG 


6 


1.2E-''= 


3 




Gene expression 




127 






10.9 


Spliceosome 


KEGG 


16 


<1E-" 




9.3 


mRNA processing 


Wikipathways 


15 


LSE-^" 


4 




Anatomical structure development 




480 






5.9 


ECM-receptor interaction 


KEGG 


15 


7.1 E-" 




3.4 


Focal adhesion 


KEGG 


25 


1.1 E-" 




6.0 


Striated Muscle Contraction 


Wikipathways 


12 


7.6E-'"* 




2.9 


Tight junction 


KEGG 


14 


1.9E-°' 




2.3 


Regulation of actin cytoskeleton 


KEGG 


17 


3.3E-''^ 




2.9 


Cell adhesion molecules (CAMs) 


KEGG 


10 


6.aE-°^ 


5 




Cell communication/immune response 




123 






13.8 


Toll-like receptor signaling pathway 


KEGG 


11 


2.0E-°' 




39 


Complement Activation, Classical Pathway 


Wikipathways 


5 


1.9E-°' 




5.0 


Chemokine signaling pathway 


KEGG 


7 


7.0E-'"' 




5.7 


Cell adhesion molecules {CAMs} 


KEGG 


5 


1.9E-°^ 




4.9 


TNF alpha Signaling Pathway 


Wikipathways 


5 


4.5E-°2 



'Differentially expressed genes were clustered using GO BP terms semantic similarity between genes as distance, to group functionally similar genes together. 

^Cluster enrichment score (ES). 

^ng, number of genes In the category. 

''Benjamini and Hochberg adjusted P value. 

doi:l 0.1 371 /journal.pone.0096491 .t004 



34]. On the other hand, CycHn D [CCm2, CCMD3) major 
regulatory proteins of the Gl phase and CDKNIB genes were 
overrepresented in SM. During Gl phase cyclin-dependent kinase 
inhibitor IB binds to cyclin D-CDK4 complexes, and thus 
promote the cell cycle arrest at Gl [35]. 

Another step of myogenesis relies on myoblast determination 
and differentiation. Determination and dififerentiation of myoblast 
progenitor is governed by the expression of a family of four muscle 
specific transcription factors called myogenic regulatory factors 
(MRFs): MYODl, MYFG, MYFS (myogenic factor 5) and MYOG 
(myogenrn) [36] . MYODl and MYF6 were specifically overrepre- 
sented in SM. During satellite cells proliferation phase, MYODl is 
expressed by activated satellite cells and governs myoblast lineage 
determination [36-38]. MYF6 is expressed in undifferentiated 
proliferating cells as well as in differentiated myoblast. Myogenic 
factor 6 seems to be implicated in both roles: myogenic 
specification genes acting on activated satellite cells and terminal 
differentiation genes [27,36,37]. 



Collectively, this set of differentially expressed genes strongly 
suggests notable diflferences in myogenic progenitor recruitment, 
proliferation and differentiation between LM and SM. LM seems 
to limit satellite cells activation through myostatin pathway. 
However, once initiated, cell cycle seems to be completed in LM 
with overrepresentation of genes involved in the control of cell 
cycle key step (Gl/S transition, S phase, S/G2 and G2/M 
transition, M phase). On the other hand, our results suggest that 
SM activates myogenic progenitors through insulin-like growth 
factor I, and that cells withdraw from cell cycle after mitosis at G 1 
phase allowing cell commitment to the differentiation program 
through MRFs expression. 

Migration is a crucial step in myogenesis as it allows myoblasts 
ahgnment before fusion in myotubes. Myoblasts specifically fuse to 
each other to form myotube, and in a second phase fuse to existing 
myotubes for muscular development, muscular maintenance or 
regeneration process [39]. The largest functional cluster underly- 
ing skeletal muscle discrepancy was related to anatomical structure 
and muscle development GO BP terms and "ECM-receptor 
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interaction" and "Focal adhesion" pathways. The extracellular 
matrix and cell adhesion molecule play an important role in 
myoblast mobility and fusion. CDH2, CD44 and CD 164 (CD 164 

molecule, sialomucin), three genes coding for cell-surface glyco- 
proteins, are involved in cell-cell interactions, cell adhesion and 
migration. CD44 and CD 164 molecules are two transmembrane 
proteins playing a key role in myoblast motility regulation [40,41]. 
Cadherins have been implicated in embryonic myoblast fusion, 
post-natal myogenesis or even regeneration process [42—44] . CA V3 
gene encodes an integral membrane protein, is induced during 
myoblast differentiation and has been implicated in myoblast 
fusion regulation and myotubes formation [45,46]. Last, adipo- 
nectin {ADIPOQ} well known for its implication in glucose 
metabolic regulation, have been also implicated in an autocrine/ 
paracrine signaling effects on myoblast differentiation and fusion 
[47,48]. These genes are overrepresented in SM. Our results 
suggest that, in SM, myoblast migration, alignment and fusion in 
myotubes are more active than in LM. Moreover, genes encoding 
giant sarcomeric protein such as JVEB and TTJ\f and genes 
encoding proteins of the sarcoplasmic reticulum membrane 
calcium release channel such as RYRl and TRD.N as well as genes 
encoding contractile proteins such MYH3, MYH8, MTH9 and 
MYHll, TKNI3 and TA\A^2 are differentially expressed between 
LM and SM. These genes are involved in terminal differentiation 
of myoblasts in sarcomerogenesis and in sarcomeric structure 
stabilization and maintenance [49-53]. This set of differentially 
expressed genes suggests that sarcomere assembly and mainte- 
nance processes are important process to characterize contrasted 
LM and SM expression profile. 

Altogether, regarding satellite cells activation, myoblast differ- 
entiation and fusion to form sarcomere, our results suggest higher 
myogenic activity in SM than in LM. 

Muscle Regeneration Process 

Last, functional enrichment analysis identified inflammatory 
immune response as relevant biological pathway to characterize 
LM and SM. SM overexpressed genes related to inflammatory 
response such as chemokine ligand (CCL21 and CCL23) and 
complement component or factor {CIQC, CIR, CIS, C3, C4A and 
CFB). In accordance, several studies have reported an induction of 
replication factors and cychns in early stage of prohferative phase 
following muscle injury. This induction is followed by upregulation 
of myogenic factor and cyclin dependant kinase inhibitors upon 
transition from proliferation phase to differentiation phase [54] , as 
well as o\ c'rcxprc'ssi()n of genes involved in inflammatory process, 
myogenic differentiation and ECM remodeling [55-57]. More- 
over, several genes products overexpressed in SM have been 
shown to be involved and specifically induced during muscle 
regeneration process. ADIPOQ^ is involved in the regenerative 
processes of skeletal muscle [47] whereas cadherins molecules are 
upregulated in activated satellite cells following injury [58]. 
Tenascin C and biglycan {BGM) two ECM glycoproteins are 
thought to be involved in muscle repair [59,60]. -SGJV mRNA 
expression which is low in mature myofibers, is highly upregulated 
during muscle regeneration in myoblast and newly formed/ 
regenerating myotubes and is concomitant with the expression of 
embryonic isofonii of myosin l)y these new myotubes [61]. 
Interestingly, SM re-expressed the embryonic (MTH3) and peri- 
natal (MTH8) isoforms of myosin heavy chain which could be 
associated with muscle regeneration. In fact, embryonic and peri- 
natal myosin isoforms disappear at birth and are progressively 
replaced by adult MHC [62,63] but re-expression of these 
developmental myosin isoform has already been reported during 
muscle regeneration [27,64]. Finally, insulin-hke growth factor I 



have been imphcated in muscle regeneration process and acting 
through a paracrine/autocrine regulation [27]. Thus, SM 
expression profile strongly suggests a regenerative muscular 
process which is characterized by expression of genes related to 
inflammatory response, fetal myogenic program and ECM 
protidns. While SM is used for locomotion, it could be therefore 
physiologically more active and subject to more minor lesions than 
the LM which is required for postural purpose. Thus, in SM, 
satellite cells might be activated and process to proliferation, 
myoblast differentiation and fusion for either muscle homeostasis 
or to form new multinucleated myotubes [27,65—67]. Concomi- 
tantiy, the mechanism of sarcomere maintenance has to incorpo- 
rate newly synthesized contractile proteins [68]. 

In conclusion, our study aimed to identify the biological events 
that underlie the differences between LM and SM metabolic and 
contractile properties by comparing their gene expression profiles. 
Results shed light on differentially expressed genes mainly related 
to myogenesis processes which suggests dissimilar post-natal 
myogenic activity between the two muscles. However we cannot 
presume if this results from dissimilar muscle maturity and/or 
from regeneration process occurring mainly in SM. This 
variability could affect muscle development and hence meat 
quality traits [13,69], thus skeletal muscle specificity should be 
taken into account to determine important features for meat 
quEility traits and identify useful biomarkers of pork quality. 

Methods 

Ethics Statement 

All samples analyzed in this study wen; collected post-mortem, 
from pigs raised and slaughtered in the context of pig meat 
production. These animals and the scientific investigations 
described herein arc therefore not to be considered as experi- 
mental animals per sc, as defined in EU directive 2010/63 and 
subsequent national application texts. Consequently, we did not 
seek ethical review and approval of this study as regarding the use 
of experimental animals. AH animals were reared and slaughtered 
in compliance with national regulations pertaining to livestock 
production and according to procedures approved by the French 
Veterinary Services. Pigs were raised on the France Hybrides 
nucleus herd of Sichamps and slaughtered on the ORLEANS 
Viandes commercial EU approved slaughterhouse according to 
standard procedures (ORLEANS Viandes, Fleury-les-Aubrais, 
France). 

Animals and Study Design 

Analyses presented here were performed on a subset of larger 
cohort, and defined as two half sibs family of 41 and 49 animals. 
The 90 pigs used in this study were non modified domestic pigs 
produced as an intercross on two successive generations between 
two terminal sire lines (FH016, Pietrain type line, and FH019 
synthetic line from Duroc, Large White and Hampshire founders, 
FRANCE HYBRIDES, St Jean de Braye, France). See Cherel 
et al. [70] for details. AH animals were raised on the same farm 
and slaughtered at an average body weight of 108 kg in the same 
commercial slaughterhouse according to standard procedures for 
commercial slaughtering (ORLEANS Viandes, Fleury-les-Aubrais, 
France). Average age was 151 days at slaughtering. Slaughtering 
method entailed exsanguination following electric stunning. The 
Longissimus and Semimembranosus muscles were sampled 20 
minutes post-mortem at the same time point following exsangui- 
nation from the same carcasses. To minimize biopsy site variation, 
each sample type (i.e. LM or SM) was collected by a single 
operator. Muscle samples were collected using a manual trocar 
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instrument and are localize in the superficial regions of the muscle. 
Muscle sample were immediately fi'ozen in liquid nitrogen. AH 
animals were genotyped as homozygous wild type genotypes NN 
and rn+rn+ with regard to the HAL and RN loci, respectively 
[71,72]. 

Total RNA Extraction 

Total RNA was extracted by crushing the fi-ozen tissue in Trizol 
reagent (Invitrogen, Cerg^'-Pontoise, France) and purified using 
Nucleospin RNA II Kit (Macherey-Nagel, Lyon, France). Total 
RNA was quantified using a NanoDrop ND-1000 spectropho- 
tometer (Thermo Scientific, lUkirch, France) and the integrity was 
assessed using the Agilent RNA 6000 Nano kit with an Agilent 
2100 Bioanalyzer (Agilent Technologies France, Massy, France). 
The RNA integrity number (RIN) was above eight for all samples. 

Microarray Design 

The "GenmascqChip", a custom 15 k pig skeletal muscle 
microarray was used in this study [15]. Microarray annotation was 
produces using BLAST 2.2.23+ [73] for megaBLAST analysis of 
the 15198 oligonucleotides scqucnci^s (60 mers) printed on the 
microarray against ENSEMBL cDNA and NCBI refSeq mammal 
databases. Annotation was based on similarity and quality criteria 
[74]. Among the 15198 probes of the GenmascqChip, 12939 
probes (i.e. 85% of the oligonucleotides) have been linked to a 
unicjuc aimotated sequence and to 9169 unique genes (i.e., 30% of 
redundancy). An 8x15 K oligo-microarray Agilent format was 
chosen, therefore one probe per microarray and eight microarrays 
were fitted in each slide. Description of the GenmascqChip is 
publicly available into the GEO repository (http:/ /www.ncbi.nlm. 
nih.gov/geo/) through GEO platform accession no. GPL11016. 

Microarray Hybridization 

Total RNA (350 ng) from each sample (90 LM, 90 SM) was 
individually labeled with Cy3 using the Low RNA Input Linear 
Amplification Kit PLUS, One-Color (ref 5188-5339, Agilent 
Technologies, Massy, France) following the manufacturer's 
instructions. Microarray hybridizations were carried out at 65°C 
for 17 hours in Agilent's SureHyb Hybridization Chambers 
containing 600 ng of Cy3-labeled cRNA sample. Slides were 
disassembled and washed in Gene Expression Wash Bufli;r 1 for 1 
minute at room temperature and then in Gene Expression Wash 
Buffer 2 for 1 minute at 37°C. Microarrays were scanned at 
5 (Xm/ pixel resolution using the Agilent DNA Microarray Scanner 
G2505B, and images were analyzed with Agilent Feature 
Extraction Software (Version 9.5), using the GE2-v5_95_Feb07 
FE extraction protocol. These MIAME compliant microarray data 
have been deposited into the GEO repository (http:/ /www.ncbi. 
nlm.nih.gov/ geo/) and are publicly available through GEO Series 
accession no. GSE33957. 

Microarray Data Analysis 

All statistical analyses were performed using R software version 
2.8.1 [75]. Raw spots intensities were first submitted to quahty 
filtration based on intensity, uniformity and saturation criteria. AH 

probes with more than 50% quality flagged spots within muscle 
were deleted from this study whereas remaining probes were 
considered significantly expressed in skeletal muscle. Processed 
signal intensities from filtered probes were natural log transformed 
and centered within sample by subtraction of the sample median 
value. Within probes, all spots that deviated by more than three 
times the standard deviation from the mean were considered as 
outliers and deleted from further analysis. AU probes whose spots 



were flagged or detected as outliers within more than 50% of 
samples per muscles were also removed from the analysis. To 
increase the robustness of differential expression analysis, probes 
with the smallest expression variability across samples were filtered 
out using K-means algorithm (k = 3) [76]. For each remaining 
probes, raw expression data were analyzed according to the 
following linear model of variance: 

Y = Sex-I-Slaughter Batch-I-Hybridization Batch within Sire-I- 
Carcass Weight+Muscle+E. 

Where Y is the raw expression data; Sex is the fixed effect of sex 
(2 levels), Slaughter Batch (2 levels) represent the effect of the 
slaughter season (summer or winter); Hybridization Batch within 
Sire represent the effect of the hybridization batch per Sire family 
(9 levels for Sire Family level 1, 11 levels for Sire Family level 2); 
Carcass Weight is a covariable to take into account the animal 
weight at slaughter time; Muscle is the fixed effect of the Muscle (2 
levels) and E is the residual. All effects exceeding the significant 
level of P<0.2, were kept in the model. The muscle effect was kept 
in each model. P-values were adjusted according to Benjamini and 
Hochberg multiple testing correction procedure [77]. Differen- 
tially expressed probes were selected using adjusted P-values of 
0.05 or less. Adjusted means for LM and SM were calculated using 
Ismeans function for the muscle effect (package Ismeans). 
Differentially expressed probes were assigned as overrepresented 
in LM or overrepresented in SM according to the greatest mean. 
Fold change value is expressed as ratio of the greatest to the least 
mean: the expression ratio of LM to SM muscle effect when genes 
are overrepresented in LM and as the expression ratio of SM to 
LM muscle effect when genes are overrepresented in SM. 

Validation of Microarray Data: Reverse Transcription and 
Quantitative Real Time PCR (RT-qPCR) 

RT-qPCR was performed using SYBR Green methodology to 
validate seven differentially expressed genes: ADAM metaUopep- 
tidase with thrombospondin type 1 motif 8 {ADAMTS8); aldolase A 
fructose-bisphosphate {ALDOA); CCAAT/ enhancer binding pro- 
tein alpha (CEBPA); muscle carnitine palmitoyltransferase IB 
[CPTlBj; diacylgly[:erol O-acyltransferase 2 {DGAT2); ryanodine 
receptor 1 (RYRl) and transforming growth factor beta 1 (TGFBT). 
Eight animals identified as non-oudiers in microarray analysis 
were randomly chosen to validate those genes. Complementary 
DNA was synthesized from 2 Hg of total RNA previously used in 
microarray analysis, using High Capacity cDNA Reverse Tran- 
scription Kit (Applied Biosystems, Foster City, CA). Primers 
(Table 5) were designed from porcine sequences using Primer 
Express software 3.0 (Applied Biosystems). Amplification was 
performed in triphcate, in 12.5 |a,l with 5 ng of reverse-transcribed 
RNA and both forward and reverse primers (200 nM each) in IX 
PCR buffer (Fast SYBR Green Master Mix, AppUed Biosystems) 
containing Uracil DNA glycosylase to prevent any DNA 
contamination from previous PCR. A StepOnePlus™ Real Time 
PCR system (Applied Biosystems) was used. Thermal cycling 
conditions were as follows: 50°C for 2 min, 95°C for 20 s, followed 
by 40 cycles of denaturation at 95°C for 3 s and annealing at 60°C 
for 30 s. Specificity of the amplification products was checked by 
dissociation curves analysis. Three genes were selected as stable 
reference genes for normalization using geNorm algorithm [78]: 
beta 2 microglobulin (B2M), TATAbox binding protein (TBP) and 
18S {18S rRNA predeveloped TaqMan kit from Applied 
Biosystems). For each sample, a normalization factor (NF) was 
calculated using geNorm algorithm and used for subsequent 
normalization. The normalized expression level (Nexp) was 
calculated according to the following formula: Ne,jp = E~^'^' 
(sample-calibrator)/ j^p^ where the Calibrator is a pool of the 16 
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Table 5. Primers pairs used in quantitative real-time RT-PCR. 





Target gene^ 


GenBank Accession^ 


Primers pair sequences (Forward/Reverse) 


ADAMTS8 


BI360058 


ACCCCTCCAGCTATGGCTACA 


TGGATGGCTCCGCTGTTT 


ALDOA 


CB286787 


CGCTGTCCCTGGGATCAC 


GCACTTGTTGATGGCGTTGA 


CEBPA 


AF1 03944.1 


GTGGACAAGAACAGCAACGA 


CTCCAGCACCTTCTGTTGAG 


CPT1B 


AF284832 


CACTGTCTGGGCAAACCAAA 


GCCACCTGGTAGGAACTCTCAAT 


DGAT2 


DT325702 


CCTGATGTCTGGAGGCATCTG 


CACGATGATGATGGCATTGC 


RYR1 


M91451 


CCCTGTGTGTGTGCAATGG 


GTTTGTCTGCAGCAGAAGCT 


TGFB1 


AF461808 


AGCGGCAACCAAATCTATGATAA 


CGACGTGTTGAACAGCATATATAAGC 




DQ178123 


AAACGGAAAGCCAAATTACC 


ATCCACAGCGTTAGGAGTGA 


TBP^ 


DQ845178 


AACAGTTCAGTAGTTATGAGCCAG 


AGATGTTCTCAAACGCTTCG 



'ADAMTS8, ADAM metallopeptidase with thrombospondin type 1 motif, 8; ALDOA, aldolase A, fructose-bisphosphate; B2M, beta-2-microglobulin; CEBPA, CCAAT/ 
enhancer binding protein (C/EBP), alpha; CPT1B, carnitine palmitoyltransferase IB (muscle); DGAT2, diacylglycerol 0-acyltransferase 2; RYR1, ryanodine receptor 1 
(skeletal); TBP1, TATA box binding protein; TGFB1, transforming growth factor, beta 1. 
^Accession number of the Sus scrofa sequence used to design primers. 
^Gene used as reference for normalization. 
doi:1 0.1 371/journal.pone.0096491 .t005 

skeletal muscle samples and E is the PGR efficiency calculated 
from the slope of calibration curve. Normalized expression levels 
of mRNAs were then compared between muscle using the Student 
t-test and P-value ^0.05 for significance. 

Functional Analysis 

To facilitate functional categorization of differentially expressed 
genes, we used hierarchical clustering of genes based on semantic 
similarity of GO BP terms. Semantic similarity was calculated 
between each pairwise combination of difFerentially expressed 
genes with a muscle fold change ratio (i.e. ratio of the greatest to 
the least muscle effect) above 1.5. Similarity was computed 
according to semantic similarity measures based on the method of 
Wang [26] implemented in the GOSemSim package [79]. Two 
genes that shared several GO BP terms result in similarity value 
close to one, indicating that they are similar in terms of biological 
process. Conversely dissimilar genes result in similarity value close 
to zero. Then hierarchical clustering was performed using "1- 
semantic similarity" as distance between two genes (similar genes 
have a distance close to zero) and "ward" as aggregation criterion 
to identify clusters of genes that share BP terms. 

Functional characterization of clustered genes was performed 
using Gene Set Analysis Toolkit V2 (WebGestalt, http;//bioinfo. 
vanderbilt.edu/webgestalt/) [80,81] using BP GO terms, KEGG 
pathways [82] and WikiPathways [83]. The lists of genes were 
uploaded using orthologous human ENTREZ gene ID. A 
minimum of five genes was required for a term to be considered 
of interest. For each terms of interest, significance levels were 
calculated following a hypergeometrical test using GenmascqChip 
15K orthologous human ENTREZ gene ID as background. A 
multiple testing correction P-value was calculated according to 



Benjamini and Hochberg procedure and an adjusted P-value of 
0.05 or less was retained for significance. 

Supporting Information 

Table SI Genes overexpressed in Longissimus. Results were 
expressed as the Longissimus to Semimembranosus ratio of the gene 
expression. The P-value of each gene was adjusted according to 
the Benjamini and Hochberg procedure. Differentially expressed 
probes were selected using adjusted P-values of 0.05 or less. 
Redundancy represented the number of probes per gene. In this 
list, 240 genes had more than one probe. 
(XLSX) 

Table S2 Genes overexpressed in Semimembranosus. Results were 
expressed as the Longissimus to Semimembranosus ratio of the gene 
expression. The P-value of each gene was adjusted according to 
the Benjamini and Hochberg procedure. Differentially expressed 
probes were selected using adjusted P-values of 0.05 or less. 
Redundancy represented the number of probes per gene. In this 
list, 413 genes had more than one probe. 
(XLSX) 

Table S3 Functional cluster composition. Annotated difFereii- 
tially expressed genes with a muscle fold change above 1.5 were 
clustered based on their functional annotation (GO BP) semantic 
similarity. Hierarchical clustering was performed using "1- 
semantic similarity" as distance between two genes (similar genes 
have a distance close to zero) to identify clusters of genes sharing 
BP terms. 
(XLSX) 
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Table S4 Enricht- d biological process of clustered differentially 
expressed genes. Functional characterization of clustered genes 
was performed using Gene Set Analysis Toolkit V2 (WebGestalt, 
http://bioinfo.vanderbilt.edu/webgestalt/) using BP GO terms. 
The lists of genes were uploaded using orthologous human 
ENTREZ gene ID. A minimum of five genes was required for a 
term to be considered of interest. For each terms of interest, 
significance levels were calculated foUowing a hypergeometrical 
test using GenmascqChip 15 K orthologous human ENTREZ 
gene ID as background. A multiple testing correction P-value was 
calculated according to Benjamini and Hochberg procedure and 
an adjusted P-value of 0.05 or less was retained for significance. 
PCLSX) 

Table S5 Enriched pathways of clustered differentially expressed 
genes. Pathway analysis of clustered genes was performed using 
Gene Set Analysis Toolkit V2 (WebGestalt, http://bioinfo. 
vanderbilt.edu/webgestalt/) using KEGG pathways and Wiki- 
Pathways. The lists of genes were uploaded using orthologous 
human ENTREZ gene ID. A minimum of five genes was required 
for a term to be considered of interest. For each terms of interest. 
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